#include <array>
#include <map>
#include <unordered_map>

#include "fit_pro.hh"
#include "TFolder.h"
#include "TFile.h"

using namespace std;

fit_pro::fit_pro(string input_filename, string align_filename, string range_filename, int layer_num, float rmse_precision, bool idx_discontinuous, int idx_disc, double sigma, double LSB, string save_path)
{
    this->det_layer_used = layer_num;
    this->rmse_precision = rmse_precision;
    this->idx_disc = idx_disc;
    this->sigma = sigma;
    this->LSB = LSB;
    this->save_path = save_path;

    read_dec_file(align_filename, alignment);
    read_dec_file(range_filename, dec_range);

    vector<float> x_hitrange(dec_range.begin(), dec_range.begin() + dec_range.size() / 2);
    vector<float> y_hitrange(dec_range.begin() + dec_range.size() / 2, dec_range.end());
    fit_range = x_hitrange;
    read_alignment();

    int data_number = read_root_file(input_filename);
    is_file_get = data_number > 0;

    setlayer = 0;
    setlayer_used = det_layer_used;
    out_result();

    idx_gen(); // prepare layer indexes for fit
    layer_len = layer_idx.size();
    convert_position_data();
    run();
}

fit_pro::~fit_pro()
{
}

int fit_pro::read_root_file(string data_infilename)
{
    hit_data_file = new TFile(data_infilename.c_str());
    if (hit_data_file->IsZombie() || !hit_data_file->IsOpen())
    {
        cout << cRED << "[ERROR]: File" << data_infilename << " open failed. Please check the file" << cRESET << endl;
        return -1;
    }
    hit_data_tree = (TTree *)hit_data_file->Get("fec_hit_data");
    if (hit_data_tree == nullptr)
    {
        cout << cRED << "[ERROR]: Cannot get TTree* fec_hit_data in file " << data_infilename << cRESET << endl;
        return -1;
    }
    hit_data_tree->SetBranchAddress("trigger_id_all", &trigger_id, &b_triggera_id);
    for (int i = 0; i < det_layer_used; i++) {

        string dec_name = "dector" + to_string(i);
        hit_data_tree->SetBranchAddress(dec_name.c_str(), &detector_data[i].x_nhits, &b_detector_data[i]);
    }
    cout << cGREEN << "TTree in TFile " << data_infilename << " get!" << cRESET << endl;
    total_entried_num = hit_data_tree->GetEntries();
    if (total_entried_num != 0)
    {
        hit_data_tree->GetEntry(0);
    }
    return total_entried_num;
}

bool fit_pro::convert_position_data()
{
    if (!is_file_get)
    {
        cout << cRED << "Data in file initial failed!" << endl;
        return false;
    }
    cout << total_entried_num << endl;

    /* convert hit data to physical value */
    for (int hit_idx = 0; hit_idx < total_entried_num; hit_idx++)
    {
        hit_data_tree->GetEntry(hit_idx);

        vector<float> temp_pos_x, temp_chny;
        vector<int> temp_amp_x, temp_ampy;
        vector<int> temp_layer_x, temp_dimy;
        vector<float> temp_clu_x, temp_cluy;

        int xnum = 0;
        int ynum = 0;

        for (int i = setlayer; i < setlayer_used; i++)
        {
            if (hit_idx == 0)
            {
                layer_z.push_back(detector_data[i].z * 10); //shouble 'z' unit as cm;
            }

            if (detector_data[i].x_nhits > 0)
            {
                xnum++;
                for (int j = 0; j < detector_data[i].x_nhits; j++)
                {
                    temp_pos_x.push_back(detector_data[i].x[j] * 0.4); //mm
                    temp_amp_x.push_back(abs(detector_data[i].x_amp[j]) * LSB); // fC
                    temp_clu_x.push_back(detector_data[i].x_chn_num[j]);
                    temp_layer_x.push_back(i);
                }
            }

            if (detector_data[i].y_nhits > 0)
            {
                ynum++;
                for (int j = 0; j < detector_data[i].y_nhits; j++)
                {
                    temp_chny.push_back(detector_data[i].y[j] * 0.4);
                    temp_ampy.push_back(abs(detector_data[i].y_amp[j]) * LSB);
                    temp_cluy.push_back(detector_data[i].y_chn_num[j]);
                    temp_dimy.push_back(i);
                }
            }

        }

        // and only keep events with all layers hitted for later usage
        /* [todo] biased selection */
        if (
            (temp_clu_x.size() > 0 && xnum == (setlayer_used - setlayer)) || 
            (temp_cluy.size() > 0 && ynum == ((setlayer_used - setlayer))))
        {
            trig_id.push_back(trigger_id);
            t_position_x.push_back(temp_pos_x);
            t_amp_x.push_back(temp_amp_x);
            t_layer_x.push_back(temp_layer_x);
            t_cluster_x.push_back(temp_clu_x);

            t_position_y.push_back(temp_chny);
            t_amp_y.push_back(temp_ampy);
            t_layer_y.push_back(temp_dimy);
            t_cluster_y.push_back(temp_cluy);
        }
        else 
          {
            // [todo] no-hit statistics
          }
        
    }
    cout << cBLUE << "Data in file convert success!" << cRESET << endl;
    cout << endl;

    return true;
}

void fit_pro::run()
{
  process_data('x');
  process_data('y');

  writeSpatial();
}

void fit_pro::process_data(char axis)
{
  /* fitting data */
  if (axis == 'x') {
      t_position = t_position_x;
      t_amp = t_amp_x;
      t_layer = t_layer_x;
      t_cluster = t_cluster_x;
      layer_select = 0;
  }
  else if (axis == 'y') {
      t_position = t_position_y;
      t_amp = t_amp_y;
      t_layer = t_layer_y;
      t_cluster = t_cluster_y;
      layer_select = 1;
  }

  /* track fit */
  data_select();
  data_fit();

  /* result */
  if (axis == 'x') {
    fec_datax = fec_data_fit;
    xtarget_data = target_data;
    xtarget_data_amp = target_data_amp;
    xtarget_data_cluster = target_data_cluster;
  }
  else if (axis == 'y') {
    fec_datay = fec_data_fit;
    ytarget_data = target_data;
    ytarget_data_amp = target_data_amp;
    ytarget_data_cluster = target_data_cluster;
  }

  //
  clear_data();
}

void fit_pro::clear_data()
{
  nr_combinations.clear();

  target_data.clear();
  target_data_amp.clear();
  target_data_cluster.clear();

  fec_data_fit.clear();
  fec_data_combine.clear();
}

void fit_pro::clear_init()
{
  trig_id.clear();
  t_position.clear();
  t_amp.clear();
  t_layer.clear();
  t_cluster.clear();
  t_position_x.clear();
  t_amp_x.clear();
  t_layer_x.clear();
  t_cluster_x.clear();
  t_position_y.clear();
  t_amp_y.clear();
  t_layer_y.clear();
  t_cluster_y.clear();
  layer_z.clear();
  layer_idx.clear();
}

bool fit_pro::data_select()
{
  if (trig_id.size() == 0) {
    cout << cRED << "no valid data" << cRESET << endl;
    cout << endl;
    return false;
  }
  else {
    /**
     * Loop all events and fill in the target layer's event buffer (ignore non-fired events)
     */
    for (int i = 0; i < t_layer.size(); i++) {
      int _max_amp = 0;
      vector<double> _hit_pos;
      vector<double> _hit_amp;
      vector<double> _cluster_size;

      /* first element is trigger id, so that each entry will be self-contained */
      _hit_pos.push_back(trig_id[i]);
      _hit_amp.push_back(trig_id[i]);
      _cluster_size.push_back(trig_id[i]);

      /* loop all clusters */
      for (int j = t_layer[i].size() - 1; j >= 0; j--) {
          if (t_layer[i][j] == idx_disc) {
              _hit_pos.push_back(t_position[i][j]);
              _hit_amp.push_back(t_amp[i][j]);
              _cluster_size.push_back(t_cluster[i][j]);
              if (t_amp[i][j] > _max_amp) {
                _max_amp = t_amp[i][j];
              }

              // remove target layer from fitting data
              t_position[i].erase(t_position[i].begin() + j);
              t_amp[i].erase(t_amp[i].begin() + j);
              t_layer[i].erase(t_layer[i].begin() + j);
          }
      }

      // append target layer data
      if (_max_amp != 0) {
        target_data.push_back(_hit_pos);
        target_data_amp.push_back(_hit_amp);
        target_data_cluster.push_back(_cluster_size);
      }
    }

    /**
     * Loop all events and fill in the nr of candidate combinations for track fit for each event
     *
     * Note: if any of the remaining layers does not contain cluster, then this
     * event has no candidate track
     */
    vector<vector<int>> containers(layer_len);
    for (int i = 0; i < t_layer.size(); i++) {
      for (int j = 0; j < layer_len; j++) {
        containers[j].push_back(count(t_layer[i].begin(), t_layer[i].end(), layer_idx[j]));
      }
    }
    for (int i = 0; i < containers[0].size(); ++i) {
      int temp = 1;
      for (int j = 0; j < layer_len; ++j) {
        temp *= containers[j][i];
      }
      nr_combinations.push_back(temp);
    }

    /**
     * Loop all events and fill in all the possible hit combinations
     */
    for (int event_id = 0; event_id < t_position.size(); event_id++) {
        vector<float> temp;
        vector<float> temp_current;
        data_combine(0, event_id, containers, temp, temp_current);
        fec_data_combine.push_back(temp);
    }

    return true;
  }
}

void fit_pro::data_combine(int layer_id, int event_id, vector<vector<int>> &containers, vector<float> &temp_x, vector<float> &temp_x_current)
{
    if (layer_id == containers.size()) {
        return;
    }

    for (int i = 0; i < containers[layer_id][event_id]; i++) {
        int k = 0;
        for (int j = 0; j < layer_id; j++) {
            k = k + containers[j][event_id];
        }

        vector<float> temp_x_next;
        temp_x_next.clear();
        if (layer_id != 0) {
            temp_x.push_back(t_position[event_id][i + k]);
            temp_x_next = temp_x_current;
            temp_x_next.push_back(t_position[event_id][i + k]);
        }
        else {
            temp_x_next.push_back(t_position[event_id][i]);
            temp_x.push_back(t_position[event_id][i]);
        }
        data_combine(layer_id + 1, event_id, containers, temp_x, temp_x_next);

        if ((i < containers[layer_id][event_id] - 1) && (layer_id != 0)) {
            temp_x.insert(temp_x.end(), temp_x_current.begin(), temp_x_current.end());
        }
    }
}

void fit_pro::data_fit()
{
  // prepare fit parametor based on target layer
  if (layer_select == 0) {
    for (int i = 0; i < det_layer_used; i++) {
      if (i == idx_disc) {
        t_target_z = layer_z[i] + t_align_z;
        layer_z.erase(layer_z.begin() + i);
      }
    }
  }

  /**
   * Loop all events with valid track candidates
   */
  vector<float> y_data;
  for (int i = 0; i < nr_combinations.size(); i++) {
    // loop all candidate tracks and select the one with min rmse
    int k = 0;
    float rmse_min = rmse_precision;
    float slope = 0;
    float intercept = 0;
    float fit_position = 0;

    /* loop all candidate combinations */
    for (int j = 0; j < nr_combinations[i]; j++) {
      float rmse = 0;
      y_data.clear();

      // prepare y-axis
      for (int m = 0; m < layer_len; m++) {
        y_data.push_back(fec_data_combine[i][k + m]);
      }

      // least-square linear fit
      Eigen::MatrixXd X(layer_len, 2);
      Eigen::VectorXd Y(layer_len);
      for (int n = 0; n < layer_len; n++) {
        X(n, 0) = 1.0;
        X(n, 1) = layer_z[n] + z_alignment[n];
        // alignment correction
        if (layer_select == 0) {
          Y(n) = y_data[n] + x_alignment[n];
        }
        else {
          Y(n) = y_data[n] + y_alignment[n];
        }
      }

      // SVD decomposition
      Eigen::VectorXd beta = X.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Y);
      for (int z = 0; z < layer_len; z++) {
        float fit_data=beta(0) + beta(1) * layer_z[z];
        rmse = rmse + (y_data[z] - fit_data) * (y_data[z] - fit_data);
      }
      rmse = sqrt(rmse / layer_len);

      // keep the one with minimum rmse
      if (rmse < rmse_min) {
          rmse_min = rmse;
          slope = beta(1);
          intercept = beta(0);
          fit_position = t_target_z * slope + intercept;
      }

      // next candidate
      k = k + layer_len;
    }

    // keep the fit result with minimum rmse, only for events with at least one track fit rmse < rmse_precision
    if (fit_position != 0 && rmse_min < rmse_precision && fit_position >= fit_range[0] && fit_position <= fit_range[1]) {
      fec_data_fit.push_back({trig_id[i], slope, intercept, rmse_min, fit_position});
    }
  }

  cout << cBLUE << "fec data fit success" << cRESET << endl;
  cout << endl;
}

void fit_pro::findTarget(
    vector<vector<double>> &aimData
    , vector<vector<double>> &ampData
    , vector<vector<double>> &clusterData
    , vector<vector<float>> &fecData
    , vector<vector<double>> &dataSecondary
    , vector<vector<double>> &dataAmpSecondary
    , vector<vector<double>> &dataClusterSecondary
    , double alignment)
{
    for (auto &row : aimData)
    {
        for (const auto &track : fecData)
        {
            if (row.front() == track.trigger_id)
            {
                auto expected_position = track.hit_position();
                for (int i = 1; i < row.size(); i++)
                {
                    row[i] = row[i] - expected_position + alignment;
                }
                dataSecondary.push_back(row);
            }
        }
    }

    for (auto &row : ampData)
    {
        for (const auto &track : fecData)
        {
            if (row.front() == track.front())
            {
                dataAmpSecondary.push_back(row);
            }
        }
    }

    for (auto &row : clusterData)
    {
        for (const auto &track : fecData)
        {
            if (row.front() == track.front())
            {
                dataClusterSecondary.push_back(row);
            }
        }
    }
}

void fit_pro::calculateSpatial(
    vector<vector<double>> &dataSecondary
    , vector<vector<double>> &dataAmpSecondary
    , vector<vector<double>> &dataClusterSecondary
    , vector<vector<double>> &dataTertiary
    , vector<vector<double>> &dataAmpTertiary
    , vector<vector<double>> &dataClusterTertiary
    , const string &filename, bool xory)
{
    for (int i = 0; i < dataSecondary.size(); i++)
    {
        vector<double> temp, amp, cluster;
        double min = sigma, ampTemp = 0, clusterTemp = 0;
        temp.push_back(dataSecondary[i][0]);
        amp.push_back(dataAmpSecondary[i][0]);
        cluster.push_back(dataClusterSecondary[i][0]);
        for (int j = 1; j < dataSecondary[i].size(); j++)
        {
            if (abs(dataSecondary[i][j]) < abs(min))
            {
                min = dataSecondary[i][j];
                ampTemp = dataAmpSecondary[i][j];
                clusterTemp = dataClusterSecondary[i][j];
            }
        }
        temp.push_back(min);
        amp.push_back(ampTemp);
        cluster.push_back(clusterTemp);
        if (temp.size() > 1 && min != sigma)
        {
            dataTertiary.push_back(temp);
            dataAmpTertiary.push_back(amp);
            dataClusterTertiary.push_back(cluster);
        }
    }

    {
        draw_spatial(dataTertiary, xory);
        draw_energy(dataAmpTertiary, xory);
        draw_cluster(dataClusterTertiary, xory);

        string eff_name;
        double eff;
        if (xory)
        {
            eff_name = "x_eff_" + to_string(idx_disc + 1);
            eff = dataTertiary.size() / total_stat_x;
        }
        else
        {
            eff_name = "y_eff_" + to_string(idx_disc + 1);
            eff = dataTertiary.size() / total_stat_y;
        }
        cout << cMAGENTA << eff_name << " " << eff << cRESET << endl;
        outall << eff_name << " " << eff << endl;
    }
}

void fit_pro::cal_all_energy(vector<vector<double>> &xampData, vector<vector<double>> &yampData, vector<double> &Ampall)
{
    for (auto &row : xampData)
    {
        for (const auto &row1 : yampData)
        {
            if (row.front() == row1.front())
            {
                double sum = row.back() + row1.back();
                Ampall.push_back(sum);
            }
        }
    }
}

void fit_pro::writeSpatial()
{
  tempx = fec_datax;
  tempy = fec_datay;
  total_stat_x = fec_datax.size();
  total_stat_y = fec_datay.size();

  vector<vector<double>> x_residual, y_residual, x_amp, y_amp, x_cluster_size, y_cluster_size;
  findTarget(xtarget_data, xtarget_data_amp, xtarget_data_cluster, tempx, x_residual, x_amp, x_cluster_size, t_align_x);
  findTarget(ytarget_data, ytarget_data_amp, ytarget_data_cluster, tempy, y_residual, y_amp, y_cluster_size, t_align_y);

  vector<vector<double>> x_min_residual, y_min_residual, x_min_res_amp, y_min_res_amp, x_min_res_cluster_size, y_min_res_cluster_size;
  calculateSpatial(x_residual, x_amp, x_cluster_size, x_min_residual, x_min_res_amp, x_min_res_cluster_size, save_path + spatial_name + to_string(idx_disc) + "idx_x_reslut.txt", true);
  calculateSpatial(y_residual, y_amp, y_cluster_size, y_min_residual, y_min_res_amp, y_min_res_cluster_size, save_path + spatial_name + to_string(idx_disc) + "idx_y_reslut.txt", false);

  vector<double> Ampall;
  cal_all_energy(x_min_res_amp, y_min_res_amp, Ampall);
  draw_all_energy(Ampall);

  outall.close();
}

void fit_pro::idx_gen()
{
    for (int i = setlayer; i < setlayer_used; i++) {
        if (i != idx_disc) {
            layer_idx.push_back(i);
        }
    }
}

void fit_pro::read_dec_file(string filename, vector<float> &dec)
{
    ifstream dec_file(filename);
    if (!dec_file.is_open())
    {
        cout << filename << " open failed" << endl;
    }
    else
    {
        float tmp;
        while (dec_file >> tmp)
        {
            dec.push_back(tmp);
        }
    }
}

void fit_pro::read_alignment()
{
  for (int i = 0; i < alignment.size(); i++) {
    if (i < det_layer_used) {
      x_alignment.push_back(alignment[i]);
    }
    else if (i < det_layer_used * 2 && i >= det_layer_used) {
      y_alignment.push_back(alignment[i]);
    }
    else if (i < det_layer_used * 3 && i >= det_layer_used * 2) {
      z_alignment.push_back(alignment[i]);
    }
  }

  for (int i = 0; i < x_alignment.size(); i++) {
    if (i == idx_disc) {
      t_align_x = x_alignment[i];
      t_align_y = y_alignment[i];
      t_align_z = z_alignment[i];
      x_alignment.erase(x_alignment.begin() + i);
      y_alignment.erase(y_alignment.begin() + i);
      z_alignment.erase(z_alignment.begin() + i);
    }
  }
}

void fit_pro::draw_spatial(vector<vector<double>> rms, bool xory)
{
    string name = (xory ? "x" : "y");
    name = name + "_spatial_layer" + to_string(idx_disc + 1);
    TH1F *spatial = new TH1F(name.c_str(), "Histogram from rms", 100, -5, 5);

    for (const auto &rms_vec : rms)
    {
        for (size_t j = 1; j < rms_vec.size(); j++)
        {
            spatial->Fill(rms_vec[j]);
        }
    }

    TF1 *fitFcn = new TF1("fitFcn", "gaus(0) + gaus(3)", -5, 5);
    fitFcn->SetNpx(1600);

    // Set initial parameter values
    fitFcn->SetParameters(spatial->GetMaximum(), spatial->GetMean(), spatial->GetRMS(),
                          spatial->GetMaximum() / 2, spatial->GetMean(), spatial->GetRMS() * 2);

    // Fit histogram
    string xais = (xory ? "#Deltax" : "#Deltay");
    spatial->Fit("fitFcn", "R");
    spatial->GetXaxis()->SetTitle(xais.c_str());
    spatial->GetYaxis()->SetTitle("Counts");
    spatial->GetXaxis()->CenterTitle();
    spatial->GetYaxis()->CenterTitle();

    double lowerBound = -5.0;
    double upperBound = 5.0;

    while (fitFcn->GetParameter(0) < 0 || fitFcn->GetParameter(3) < 0)
    {
        lowerBound += 0.5;
        upperBound -= 0.5;
        delete fitFcn;
        fitFcn = new TF1("fitFcn", "gaus(0) + gaus(3)", lowerBound, upperBound);
        fitFcn->SetNpx(1600);
        fitFcn->SetParameters(spatial->GetMaximum(), spatial->GetMean(), spatial->GetRMS(),
                              spatial->GetMaximum() / 2, spatial->GetMean(), spatial->GetRMS() * 2);
        spatial->Fit("fitFcn", "R");
    }

    TCanvas *c2 = new TCanvas("c2", "spitial_resolution", 800, 600);

    // Draw original histogram and fit function on new canvas
    spatial->Draw();
    // Create TF1 objects representing two Gaussian functions
    TF1 *gaus1 = new TF1("gaus1", "gaus", -5, 5);
    TF1 *gaus2 = new TF1("gaus2", "gaus", -5, 5);

    // Set their parameters
    gaus1->SetParameters(fitFcn->GetParameter(0), fitFcn->GetParameter(1), fitFcn->GetParameter(2));
    gaus2->SetParameters(fitFcn->GetParameter(3), fitFcn->GetParameter(4), fitFcn->GetParameter(5));
    gaus1->SetLineColor(kBlue);
    int brown = TColor::GetColor("#8B4513");
    gaus2->SetLineColor(brown);
    // Draw them on the same canvas
    gaus1->Draw("same");
    gaus2->Draw("same");

    double sigma = (fitFcn->GetParameter(0) > fitFcn->GetParameter(3)) ? fitFcn->GetParameter(2) : fitFcn->GetParameter(5);
    cout << cMAGENTA << "Sigma of Gaussian: " << sigma << endl;
    auto t = new TLatex;
    t->DrawLatex(1.25, 30, to_string(sigma).c_str());

    // Save canvas as PDF file
    string filename;
    string sigmaname;

    if (xory)
    {
        filename = save_path + to_string(idx_disc) + "_spatial_x" + ".png";
        sigmaname = "Xspatial_" + to_string(idx_disc + 1);
    }
    else
    {
        filename = save_path + to_string(idx_disc) + "_spatial_y" + ".png";
        sigmaname = "Yspatial_" + to_string(idx_disc + 1);
    }

    outall << sigmaname << " " << sigma << endl;
    c2->SaveAs(filename.c_str());

    spatial->Reset();
    delete spatial;
    delete fitFcn;
    delete gaus1;
    delete gaus2;
    delete c2;
    delete t;
}

void fit_pro::draw_energy(vector<vector<double>> energy, bool xory)
{
    string name = (xory ? "x" : "y");
    name = name + "_energy_layer" + to_string(idx_disc + 1);

    auto* view_fout = new TFile("fit_pro.root","update");
    view_fout->cd();

    TH1F *muonenergy = new TH1F(name.c_str(), "Histogram from muon deposit energy", 100, -5, 800);

    for (const auto &energy_vec : energy) {
        for (size_t j = 1; j < energy_vec.size(); j++) {
            muonenergy->Fill(energy_vec[j]);
        }
    }

    TF1 *fitFcn = new TF1("fitFcn", "[0]*TMath::Landau(x,[1],[2])*TMath::Gaus(x,[3],[4])", -5, 800);
    fitFcn->SetNpx(500);

    // 设置初始参数值
    fitFcn->SetParameter(0, muonenergy->GetMaximum()); // 高度
    fitFcn->SetParameter(1, muonenergy->GetMean());    // 朗道的中心
    fitFcn->SetParameter(2, muonenergy->GetRMS());     // 朗道的宽度
    fitFcn->SetParameter(3, muonenergy->GetMean());    // 高斯的中心
    fitFcn->SetParameter(4, muonenergy->GetRMS());     // 高斯的宽度

    // 拟合直方图
    muonenergy->Fit("fitFcn", "RQ");


    muonenergy->Write();

    // 创建一个新的画布
    TCanvas *c2 = new TCanvas("c2", "Canvas", 800, 600);

    // 在新的画布上绘制原始直方图和拟合函数
    muonenergy->Draw();
    fitFcn->Draw("same"); // 在原始直方图上添加拟合函数的图像

    double energy_mip = fitFcn->GetParameter(1);
    auto t = new TLatex;
    t->DrawLatex(500, 200, to_string(energy_mip).c_str());
    cout << cMAGENTA << "mip of land-Gaussian: " << energy_mip << endl;

    // 保存画布为PDF文件
    // Save canvas as PDF file
    string filename;
    string energyname;

    if (xory)
    {
        filename = save_path + to_string(idx_disc) + "_energy_x" + ".png";
        energyname = "Xmuonenergy_" + to_string(idx_disc + 1);
    }
    else
    {
        filename = save_path + to_string(idx_disc) + "_energy_y" + ".png";
        energyname = "Ymuonenergy_" + to_string(idx_disc + 1);
    }

    outall << energyname << " " << energy_mip << endl;
    c2->SaveAs(filename.c_str());

    muonenergy->Reset();
    delete muonenergy;
    delete fitFcn;
    delete c2;


    view_fout->Write();
    view_fout->Close();
}

void fit_pro::draw_all_energy(vector<double> energy)
{
    string name = "xy_energy_layer" + to_string(idx_disc + 1);
    TH1F *muonenergy = new TH1F(name.c_str(), "Histogram from muon deposit energy", 100, -5, 1200);

    for (const auto &energy_vec : energy) {
        muonenergy->Fill(energy_vec);
    }

    TF1 *fitFcn = new TF1("fitFcn", "[0]*TMath::Landau(x,[1],[2])*TMath::Gaus(x,[3],[4])", -5, 1200);
    fitFcn->SetNpx(500);

    // 设置初始参数值
    fitFcn->SetParameter(0, muonenergy->GetMaximum()); // 高度
    fitFcn->SetParameter(1, muonenergy->GetMean());    // 朗道的中心
    fitFcn->SetParameter(2, muonenergy->GetRMS());     // 朗道的宽度
    fitFcn->SetParameter(3, muonenergy->GetMean());    // 高斯的中心
    fitFcn->SetParameter(4, muonenergy->GetRMS());     // 高斯的宽度

    // 拟合直方图
    muonenergy->Fit("fitFcn", "R");

    // 创建一个新的画布
    TCanvas *c2 = new TCanvas("c2", "Canvas", 800, 600);

    // 在新的画布上绘制原始直方图和拟合函数
    muonenergy->Draw();
    fitFcn->Draw("same"); // 在原始直方图上添加拟合函数的图像

    double energy_mip = fitFcn->GetParameter(1);
    auto t = new TLatex;
    t->DrawLatex(600, 100, to_string(energy_mip).c_str());
    cout << cMAGENTA << "XYmip of land-Gaussian: " << energy_mip << endl;

    // 保存画布为PDF文件
    string filename = save_path + to_string(idx_disc) + "_energy_xy.png";
    string energyname = "XYmuonenergy_" + to_string(idx_disc + 1);
    outall << energyname << " " << energy_mip << endl;
    c2->SaveAs(filename.c_str());

    muonenergy->Reset();
    delete muonenergy;
    delete fitFcn;
    delete c2;
}

void fit_pro::draw_cluster(vector<vector<double>> cluster, bool xory)
{
    string name = (xory ? "x" : "y");
    name = name + "_hit_strip_num_layer" + to_string(idx_disc + 1);
    TH1F *stripnum = new TH1F(name.c_str(), "Histogram from muon hit strip  number", 30, 0, 30);

    for (const auto &cluster_vec : cluster)
    {
        for (size_t j = 1; j < cluster_vec.size(); j++)
        {
            stripnum->Fill(cluster_vec[j]);
        }
    }

    // Create a new canvas
    TCanvas *c1 = new TCanvas("c1", "Canvas", 800, 600);

    // Draw the histogram on the canvas
    stripnum->Draw();

    // Save the canvas as a PNG file
    string filename;
    if (xory)
    {
        filename = save_path + to_string(idx_disc) + "_stripnum_x.png";
    }
    else
    {
        filename = save_path + to_string(idx_disc) + "_stripnum_y.png";
    }
    c1->SaveAs(filename.c_str());

    // Clean up
    stripnum->Reset();
    delete stripnum;
    delete c1;
}

void fit_pro::out_result()
{
    string name = "xyunsim_";

    outall.open(save_path + name + "result.txt", ios::app);
    if (!outall.is_open())
    {
        cout << cRED << "Failed to open file: " << (save_path + name + "result.txt") << endl;
    }
}
